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The nature of dark matter is unknown. A number of dark matter candidates are quantum flavor- 
mixed particles but this property has never been accounted for in cosmology. Here we explore this 
possibility from the first principles via extensive A^-body cosmological simulations and demonstrate 
that the two-component dark matter model agrees with observational data at all scales. Substantial 
reduction of substructure and flattening of density profiles in the centers of dark matter halos found 
in simulations simultaneously resolve several outstanding puzzles of modern cosmology. Predictions 
for direct and indirect detection dark matter experiments are also made. 



PACS numbers: 95.35. +d, 95.30.Cq, 98.80.-k, 14.80.-j 

Dark matter (DM) constitutes about 80% of matter 
and 25% of the total energy density in the universe but 
its nature remains completely unknown. The existence 
of DM requires revision of the present day physics. Most 
likely, DM is a hypothetical particle or particles beyond 
the standard model [T]. 

The current heuristic paradigm of cold dark matter 
with a cosmological constant (ACDM) is very success- 
ful at reproducing the large-scale structure of the uni- 
verse but appears to disagree with observations on small 
scales. First, simulations predict the overabundance of 
small mass (dwarf) halos as compared to the much lower 
number of the observed satellite galaxies in the Local 
Group [IHS] and in the field as inferred from ALFALFA 
survey This problem was termed the "substructure 
problem" or the "missing satellite problem." Second, 
the cuspy p cx r^-^ DM density profiles in ACDM sim- 
ulations ff disagree with rotation curves of dwarf and 
low surface brightness galaxies, which indicate flattened 
or cored density profiles [5HTU]. Observations of galaxy 
clusters also indicate the presence of cores [TT]. More- 
over, the largest CDM subhalos in the Local Group-type 
environments are too dense in their centers to host any 
of the dwarf spheroidal galaxies around the Milky Way 
or Andromeda galaxies [HI [13]. These two, perhaps re- 
lated, "inner halo problems" are known as the "core-cusp 
problem" and the "to big to fail problem" , respectively. 
Importantly, numerous attempts to reconcile the ACDM 
model with observations using baryonic processes made 
so far (modified star formation, tidal gas stripping, su- 
pernova feedback) are unconvincing [H [T^HTB] . This is 
because the inner halo problems require strong feedback 
and, hence, larger star formation whereas the missing 
satellite problem requires the star formation to be sup- 
pressed. Contrary to early expectations, a mild modifica- 
tion of the DM paradigm to adopt the warm DM model 
[TOl also fails to resolve all these problems altogether 
[^T| - f25] due to similar constraints on the DM mass. 

The inability of conventional physics to resolve 
aforementioned problems within the coUisionless CDM 
paradigm indicates that cosmological DM exhibits non- 
gravitational properties as well. The most natural alter- 



native is to admit a large interaction cross-section of DM 
with itself [24l [25] but not with normal matter. Contrary 
to early claims [26l - f28] , such self- interacting dark matter 
(SIDM) was successful to explain the origin of cores pro- 
vided certain constraints on the velocity-dependent cross- 
section a{v) are satisfied [^5H55] : however it completely 
fails to solve the missing satellite problem [33]. On the 
positive side, SIDM also naturally predicts the very rapid 
formation of supermassive black holes [3S1 [37] via gravi- 
tational collapse of central (coUisional) parts of DM ha- 
los [35| [33] — the process that is absent in the "vanilla 
CDM" paradigm. At last, the existence of the narrow 
plane of Andromeda dwarf satellites [40], which has no 
explanation within coUisionless CDM, may potentially be 
addressed in SIDM, because coUisionality induces viscous 
friction on subhalo motions in massive halos (whether its 
magnitude is enough remains to be seen). 

However, an important possibility has been entirely 
missing from consideration so far, namely that a number 
of DM candidates are quantum flavor-mixed particles, 
e.g., a neutralino, an axion, a sterile neutrino, and oth- 
ers. In this Letter we demonstrate from the first prin- 
ciples via state-of-the-art iV-body cosmological simula- 
tions that even the simplest model with two-component 
quantum-mixed DM with small mass-degeneracy agrees 
with observational data at both large and small scales, 
thus settling the above problems altogether. Moreover, 
it also agrees with observational constraints on the cross- 
section set by other coUisional DM models [3TH35] . At 
last, the model makes predictions for and is testable with 
direct and indirect detection DM experiments. 

The physical idea is very simple. A mixed particle of 
a particular flavor a is a superposition of several mass- 
eigenstates = ai |mi) -I- 02 |m2) -I- . . . , where |/) 

and I to) denote wave- functions being flavor and mass 
eigenstates, and Oi, 02,... are complex constants be- 
ing elements of a unitary matrix. The dynamics of non- 
relativistic mixed particles has been shown to be rather 
unusual in that elastic scattering of mixed particles can 
lead to their irreversible escape from the gravitational po- 
tential well [m |42] leaving nothing behind. This occurs 
because scatterings of particle mass eigenstates cause 
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their inter-conversions which change their kinetic (but 
not total) energy [41], thus giving a kick of the order of 
Vk — cy^2Am/m, where Am is the mass difference of the 
mass eigenstates and c is the speed of hght. This effect, 
referred to as the "quantum evaporation" (see Ref. [12] 
for a full quantum mechanical description) and suggested 
to occur in DM halos, can simultaneously soften the den- 
sity cusps and reduce the number of subhalos. However, 
the structure on large scales where the escape velocities 
are much larger than is unaffected. 

Here we consider the simplest model of a flavor-mixed 
DM particle involving only two mass eigenstates and, cor- 
respondingly, two flavors [H] |42] , i.e. the two-component 
dark matter (2cDM) model. The masses of the mass 
eigenstates are mh and m/ < to/j, referred to as 'heavy' 
and 'light'. Since mass eigenstates have different veloci- 
ties, they propagate along different geodesies so they can 
be spatially separated by gravity during structure forma- 
tion, when the characteristic (escape) velocities are of the 
order oivk- Thereafter, DM halos are self-gravitating en- 
sembles of non-overlapping heavy and light eigenstates. 

A collisions of mass eigenstates can lead to either scat- 
tering or conversion, such as Inih} — >■ {mi), or simply 
h ^ I. Important conversions are hh — >■ hi, hh — >■ II and 
hi — ^ II (in which one or two heavy states are converted 
into the lighter states) because they increase particles' 
velocities and lead to the quantum evaporation. Because 
of energy-momentum conservation, the kinetic energy of 
the eigenstates in heavy-to-light conversions increases by 
Amc^ = (m/i — TO/)c^ in processes like hh — )■ hi and twice 
as much in hh — )• U. The reverse processes hi — > hh, 
11 — hi and 11 — > hh and can also occur if kinemati- 
cally allowed, i.e., if kinetic energy is large enough to 
produce a heavy eigenstate. Finally, the scattering pro- 
cesses II — !> II, hi — ^ hi and hh — ?> hh can occur as well, 
but there is always a set of parameters (the interaction 
strengths and the mixing angle) when they are vanishing 
[4T] |42] , hence they can be omitted in simulations when 
needed. Generally, any velocity-dependent cross-section, 
<j(v), consistent with SIDM constraints, can equally well 
be utilized in the 2cDM model. Here we use the simplest 
one in which the conversion cross-section at low veloci- 
ties is a{v) cx l/w" with a = 1 (or greater) [42]. Such 
l/i;-dependence is not ruled out by observational data 
[551 ISIl 132] , so it is used in simulations reported here and 
a/m = 0.75 cm^/g at v — Vk (for all conversion types) is 
also within the allowed range 0.1 < a/m < 0(1) cm^/g 
at w ^ 100 km/s [29l |3TJ |39] , where m is the characteristic 
mass of a DM particle, m ~ m/j. 

The physics of mixed-particle interactions was imple- 
mented in the cosmological TreePM/SPH publicly avail- 
able code [35 Gadget-2. We simulate two types of DM 
particles representing h and / mass eigenstates; the total 
numbers of each type can change due to particle con- 
versions. In order to implement interactions of DM par- 
ticles, they are treated in the code as smooth-particle- 



hydro (SPH) particles but without hydro-force accelera- 
tion. To model particles' binary interactions, we use the 
Monte-Carlo technique together with the "binary colli- 
sion approximation" [5S1 [55] , which is reliable for weakly 
coUisional systems. The Monte-Carlo particle interaction 
module works as follows. For each randomly chosen pro- 
jectile particle, Si, a nearest neighbor is found; this is 
the target particle, U (the subscript i stands for 'initial', 
i.e. before the interaction). This procedure identifies the 
type and velocity of each particle in the interacting pair, 
i.e., the input channel. For each input channel, Siti there 
always are four output channels, Sftf (here / stands for 
'final', i.e., after the interaction), namely: hh, hi, Ih 
and II. We compute the probabilities for each process 
Siti — >■ Sftf as follows: 

Ps,u^s;tf = {ptjmt,) as,ti~^s}ti\vti - v^JAt Q{Esftj) 

(1) 

where as-t^^sftf = cr is the cross-section which is as- 
sumed to be the same for all conversion channels and 
zero for scattering channels for 2cDM and vice versa for 
SIDM, vt . — v^. is the relative velocity of particles in the 
pair, pt- is the density of target species. At is the iter- 
ation time-step and Q{Es^tf ) is the Heaviside function 
which ensures that the process is kinematically allowed 
(i.e., negative final kinetic energy, E^^tf < 0, means the 
process cannot occur). Generally, as^ti^sftf depends on 
the mixing angle 9, so the above choice of cr's corresponds 
to maximal mixing 9 = tt/A for 2cDM and no mixing 
9 = ior SIDM 42J . The densities pt- of each species at 
each particle's position are computed using the appropri- 
ately modified density routine used in the SPH module 
of the original code. Whether an interaction occurs and 
through which channel it proceeds is determined by ran- 
dom drawing in accordance with the computed probabil- 
ities. If an interaction of any kind occurs, the pair kine- 
matics is computed in the center of mass frame, where 
the momentum is conserved manifestly. If a scattering 
occurs, the particles are given random antiparallel veloc- 
ities (in the center of mass frame) with magnitudes set by 
the energy and momentum conservation laws. If a con- 
version occurs, then (i) the type of one or both particles 
is changed accordingly, (ii) the magnitudes of the final 
velocities are computed accounting for the Amc^ given 
or taken, depending on the type of conversion and (iii) 
these velocities are assigned to particles in antiparallel 
directions in the center of mass frame. If no interaction 
occurs, the particle velocities and types remain intact. 
After this, the pair is marked inactive until the next time- 
step. This process is repeated for all active particles at 
each time step. 

Simulations reported here were performed using 
XSEDE high performance computing systems Trestles 
and Ranger. The 2cDM nms have 2 x 400^ = 128 mil- 
lion SPH-DM particles (in 2cDM, the initial numbers of 
h and I particles are equal) in the box of 50/i^^ Mpc 



FIG. 1. Dark matter distribution in a region of size 5h-^ Mpc 
with the standard ACDM (left panel) and the 2cDM (right 
panel) model with a/m = 0.75 cm^/g, Vk ~ 50 km/s. Both 
these simulations have the same number of particles. Note 
the deficit of substructure (small yellow clumps) in the 2cDM 
vs ACDM model. 



(comoving) with the force resolution scale of 3.5/i~^ kpc, 
and the reference ACDM run has 2 x 640^ ~ 524 mil- 
lion particles and the force resolution of 2.2h~^ kpc. Our 
box size was optimized to be large enough to be a rep- 
resentative sample the universe volume, yet it provides 
reasonable resolution at small scales. All the runs are 
DM-only simulations using the standard cosmological pa- 
rameters ftm ~ 0.3, il\ — 0.7, fib — and h — 0.7. Initial 
conditions are generated using N-GenIC code with the 
Eisenstein-Hu spectrum model, with erg — 0.9 and the 
initial redshift z = 50. Post-processing was done with 
AHF code [H], which was used to construct the halo 
mass and velocity functions, analyze halo density pro- 
files, etc. Simulations of SIDM have also been done for 
the same cosmological parameters. They fully confirm 
earlier studies, e.g., the inability to resolve the substruc- 
ture problem, hence these results are not reported here. 
A number of simulations were performed to explore a 
range of the 2cDM model parameters Am/m and a/m, 
to compare with the reference CDM and SIDM models 
and to check for numerical convergence. All the results 
will be reported in detail elsewhere; here we show the 
most important ones. 

Simulations with large mass difference nih > mi (not 
presented here) show the enhancement of substructure 
relative to the large-scale halos (masses large halos are 
diminished to produce more smaller ones). This makes 
the substructure problem worse and, hence rules out the 
non-degenerate case. We focus on the degenerate case 
mi ~ mil — m hereafter. 

The representative case of the DM distribution at 
z = for the 2cDM and ACDM models are presented in 
Fig. [l] which show the zoomed-in region of 5 Mpc across. 
One can see the reduced number of subhalos and the less 
concentrated central parts in the 2cDM case. The 2cDM 
parameters used are Am/m ~ 10~®, which corresponds 
to Vk = 50 km/s, and a/m = 0.75 cm^/g at u = Vk, 
which is fully consistent with observational constraints 
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FIG. 2. The maximum circular velocity function for the stan- 
dard ACDM model (black curve) and the 2cDM model with 
Vk = 50 km/s and a/m = 0.75 cm^/g (blue curve). The 
2cDM provides an excellent fit to the appropriately scaled 
Local Group data [21 14] (magenta points) . 

on the SIDM cross-section [55H311I3S]- For these values, 
the 2cDM circular velocity function matches the Local 
Group data the best, as is illustrated in Fig. [2] This 
figure shows the number of halos with the maximum cir- 
cular velocity above a certain value, N{> Vc^may/) versus 
Vc,max, for 2cDM and ACDM; the data points are from 
The amount of substructure is volume-dependent, 
so we appropriately renormalized the data points to re- 
produce the results of Refs. [21 14] using the velocity func- 
tion from our ACDM simulation; the procedure is legiti- 
mate for a scale- free ergodic distribution of DM substruc- 
ture. Scanning though the 2cDM model parameters, we 
have found that Vk uniquely determines the position of 
the break in the velocity function, V^^^}^ ~ Vk, whereas 
a/m determines the slope below the break. Comparison 
with observations determines Vk rather accurately to be 
around ^ 50 — 70 km/s. Interestingly, a similar value 
of a characteristic velocity < 100 km/s was found in an- 
other independent analysis of survey data ^5^ The 'best 
fit' cross-section is a/m ^ 0.75 cm^/g at Vk but values 
a factor of two smaller or larger are marginally accept- 
able too. The halo mass function exhibits even sharper 
break at halo masses M ~ 10^" A/©. Thus, the overall 
suppression of the subhalo abundance at small velocities 
and masses resolves the "missing satellite problem" . 

Fig. [3] shows the ACDM and 2cDM halo density pro- 
files. Plotted are over one hundred of well-resolved pro- 
files, the inner parts of them were truncated according to 
the numerical binary collision criterion [45] . hence they 
are trustworthy everywhere. The ACDM profiles agree 
with the NFW profile. In contrast, the 2cDM inner pro- 
files are less centrally concentrated and shallower (and 
some are core-like), as is also seen from Fig. |4] Here, 
the effective power-law index is obtained by fitting the 
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FIG. 3. Density profiles of 120 well-resolved dark halos are 
plotted for the classical ACDM (left panel) and two 2cDM 
models. The profiles are color-coded by the halo mass: red - 
most massive, blue - less massive. 
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FIG. 4. Histograms of the slopes of the inner density profiles 
of the halos shown in Fig. [S] Whereas the CDM profiles show 
a cusp with a ~ — 0.8 . . . — 1 consistent with earlier studies, 
the 2cDM profiles are much shallower: a ~ —0.2 ... — 0.6. 



profiles with the function p = po r°'{l + r/rc)^ and then 
evaluating a at r = 7 kpc/ft,. The distribution of the 
slopes ranges within a ~ —0.8 ... — 1 for CDM indi- 
cating a cusp and within a ~ —0.2 ... — 0.6 for 2cDM, 
which thus explains the "core-cusp problem" and, likely, 
the "too big to fail problem" . However, there is large di- 
versity in the values of the central density and the inner 
slopes due to different formation histories in agreement 
with a 'counter-example' study [5T]. Importantly, density 
profiles and the core sizes of massive halos are mostly sen- 
sitive to the value oia/m, the value of Vk plays little role, 
if any. The profiles of smaller halos with masses around 
or below IO^'^Mq can be sensitive to Vk- Dedicated high- 
resolution simulations are needed to study this in detail. 

The 2cDM theory is testable with direct detection ex- 
periments. Indeed, DM is a collection of h and I species 
ant they can convert into one another when they interact 
with normal matter in a detector. These conversions will 
result in energy 'mismatch' of ~ Amc^, i.e., the events 
will look like inelastic collisions. Down-conversions h I 
should always occur, along with elastic scatterings, and 
will result in additional recoil energy. In contrast, up- 
conversions I — >■ h will result in the deficit of energy and 



can occur only if kinetic energy is enough. Since Vk is 
comparable to the velocity of Earth around the Sun and 
of the Sun in the Galaxy, annular modulation of con- 
version rates is expected. The absolute value of Am is 
not constrained by our analysis, but only the relation 
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If inelastic effects are responsible for 



DAMA and CoGeNT anomalies, then Am ~ few keV. 
Then the DM mass should be mx ^ few x 10^ GeV, 
which consistent with mx ^ 130 GeV from possible an- 
nihilation gamma-ray line found in Fermi LAT data |46) . 

To summarize, cosmological A'^-body simulations with 
two-component DM have been performed for the first 
time and show different structure formation on small 
scales. The 2cDM model does not change the linear 
power spectrum (unlike warm DM); all the changes oc- 
cur in the nonlinear stage. The results are as follows. 
First, cosmology with at least two mixed DM species 
with high mass degeneracy, Am/m ^ 10^^, is preferred. 
In turn, self-interacting DM models with single species 
(e.g., non-flavor-mixed) candidates such as gravitino and 
others and non-degenerate multi-component DM models 
are disfavored. Second, 2cDM has a reduced amount of 
dwarf halos with masses below lO^^M© and typical ve- 
locities smaller than 50-70 km/s, hence it resolves the 
"substructure" problem. Third, 2cDM softens central 
cusps and makes halos less centrally concentrated, thus 
solving the "core-cusp" and "too big to fail" problems in 
the same way SIDM does. 2cDM agrees with observa- 
tions within a range of velocity-dependent cross-section 
allowed for SIDM [29H34l[39]. The constraints are tight: 
if a is too small, then it is cosmologically uninterest- 
ing, if it is too large, then the cusps will be enhanced 
due to gravithermal collapse of halos. This fine tuning, 
rephrased as "Why now?" question is a caveat of 2cDM 
model. However, SIDM and dark energy /cosmological 
constant face the same problem. Whether they are re- 
lated remains to be seen. Fourth, being a coUisional 
model, akin to SIDM, 2cDM can also explain the ex- 
istence of supermassive black holes at high-z via semi- 
coUisional DM accretion [3S]. Fifth, 2cDM passes the 
early universe constraint - that at high densities higher 
mass particles must convert into the lightest one - be- 
cause the conversion cross-section in flat space-time is 
suppressed by (Am/m)^ ~ 10~^^ over it's current value 
|42) . Sixth, conversions can be seen in direct detection 
DM experiments as inelastic (endothermic or exothermic) 
processes with AE ~ ±Atoc^ or twice as much. Seventh, 
we predict that evaporation of dwarf halos should enrich 
the environment with metals from stars formed earlier 
in these halos, which would mimic the effect of super- 
nova feedback. Since not all the small halos are evap- 
orated, the residual substructure can be responsible for 
flux anomalies in gravitational lensing observations. 
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